2 4 Ju n 19 99 Phase Diagram of the Lattice Restricted Primitive Model

We present a comprehensive study of the lattice restricted primitive model, i.e., a lattice gas consisting of an equal number of positively and negatively charged particles interacting via on-site exclusion and a 1/r potential. On the cubic lattice, Monte Carlo simulations show a line of Néel points separating a disordered, high-temperature phase from a phase with global antiferromagnetic order. At low temperatures the (high-density) ordered phase coexists with the (low-density) disordered phase. The Néel line meets the coexistence curve at a tricritical point, Tt ≃ 0.14, ρt ≃ 0.4. A simple mean-field analysis is in qualitative agreement with simulations. INTRODUCTION It has long been realized that the presence of an appreciable concentration of free ions in a system of overall charge neutrality gives rise to a variety of thermodynamic features not found in un-ionized systems. For example, in the restricted primitive model (RPM), a system of charged hard-sphere anions and cations, all with equal charge magnitude |q| and diameter σ, one has the famous limiting law established by the work of Debye and Hückel [1] in 1923, − A ex kBTV → Γ 3 D 12π as ρ → 0, T fixed. (1) Here A is the Helmholtz free energy (excess to the ideal-gas condition), V is the system volume, kB Boltzmann’s constant, T the temperature, ρ the total number density of anions and cations, and ΓD is the inverse Debye length, ΓDσ 2 = 4πρq/kBT . This is in marked contrast to an un-ionized fluid, for which 1) electronic address: [email protected] 2) On leave of absence from: Department of Physics and Astronomy, Herbert H. Lehman College, City University of New York, Bronx, NY, 10468-1589. 3) electronic address: [email protected] A/kBTV approaches a second virial-coefficient term proportional to ρ 2 rather than ρ. Results pertaining to phase transitions and criticality in the RPM did not emerge until much after the work of Debye and Hückel. It was not until 1976 that a systematic study by Stell, Wu, and Larsen appeared with the conclusion that the RPM should exhibit a liquid-gas coexistence curve with a critical point [2]. They further concluded that the critical density is much lower than for a simple argon-like fluid model, but found that that the shape and location of the coexistence curve was very sensitive to the details of the approximations that must inevitably be used in such a theoretical study. Monte Carlo simulations have yielded results fully consistent with their conclusions, but only in the last few years have the studies of different groups converged toward common values for the critical density and temperature of the RPM [3–5]. As one of us has discussed in earlier studies [6,7], the RPM can be thought of as a spin system with a long-ranged antiferromagnetic interaction J(r) = − 2 rǫ , (2) and one-dimensional spins si with |si| = 1; φij(r) = −J(r)si · sj . (3) (Since the “spins” are charges, we should perhaps refer to the interaction as antiferroelectric rather than antiferromagnetic.) The lattice-gas version of such a spin fluid is simply the spin-1 Ising model with |s| = 1, 0, or -1, with s = 0 representing a vacant site. We shall refer to this model, which forms the subject of this report, as the lattice restricted primitive model (LRPM). For a lattice gas with a nearest-neighbor J(r) instead of a Coulombic J(r), the spin-1 Ising model is a Blume-Capel model [8,9], which becomes equivalent, at full occupancy, to a nearest-neighbor spin/2 Ising model that is exactly-soluble on many two-dimensional lattices [10,11]. On the square lattice one finds a Néel point as one lowers the temperature in the absence of an external field. Moreover, although exact results are lacking in three dimensions, one expects for bipartite lattices such as the simple cubic, that as one lowers the density along the λ-line of Néel points in the ρ-T plane, one encounters a tricritical point below which the lattice gas phase-separates into paramagnetic and antiferromagnetic phases. This immediately raises the question whether this remains true for a lattice gas with a Coulombic J(r). Several years ago, in order to better understand the LRPM, we initiated Monte Carlo simulations on a simple cubic lattice. Our results are fully consistent with the presence of a λ-line of Néel points terminating in a tricritical point, below which (in T ) two phases coexist. Preliminary results were reported in Ref. [12]; here we provide more extensive results as well as details of our simulation procedure. We have supplemented our simulations with a mean-field analysis that exploits the similarity between the LRPM and a nearest-neighbor lattice gas. Høye and Stell [13] have also made a more general study of the spin-1 Ising model that included the solution of the mean-spherical approximation (MSA) for a range-parametrized J(r). As discussed by these authors, the MSA is not appropriate to study criticality and phase separation in the spin-1 model with Coulombic J(r) (although, as they show, it does yield the correct Debye-Hückel limiting law). The Høye-Stell study, however, goes beyond the MSA with a number of general observations that strongly support tricriticality for a Coulombic J(r). More recent studies made by Ciach and Stell [14] of ionic models that include the continuum-space and lattice RPM reveal that generically such models can be expected to manifest criticality, tricriticality, or both, depending upon the precise forn of their Hamiltonians. For example, the work suggests that some extendedcore lattice models may well show both liquid-gas like criticality and, at higher densities, a λ-line with associated tricriticality. For the RPM, the critical properties that follow from application of RG analysis are found to be in the Ising universality class. MODEL We consider a lattice gas of particles interacting via site exclusion (multiple occupancy forbidden) and a Coulomb interaction u(rij) = sisj/rij , where rij = |ri − rj| is the distance separating the particles (located at lattice sites ri and rj), and si = +1 or −1 is the charge of particle i. Exactly half the N particles are positively charged; the remainder are negative. They are restricted to a simple cubic lattice of V = L sites, with periodic boundaries. We assume a lattice constant a of unity, and adopt units in which q/akB = 1, q being the magnitude of the charge. In what follows we treat charge, length, energy and temperature as dimensionless. We note in passing that in lattice simulations of ionic systems an alternative definition of the potential is possible, i.e., u(rij) = sisjΦ(rij), where Φ is the Green’s function for Poisson’s equation on the lattice in question. While the simple 1/r potential and the lattice Green’s function differ somewhat at short distances, they have the same asymptotic (large r) behavior, and we expect (qualitatively) the same phase diagram in either case. In this study we restrict attention to the simple cubic lattice. This lattice, like the body-centered cubic, admits a decomposition into two sublattices (the sites on one sublattice having all nearest neighbors in the other sublattice), which obviously facilitates formation of an ordered state resembling an ionic crystal. Indeed, it was proven some time ago that the LRPM on the simple cubic lattice exhibits longrange order at sufficiently low temperatures and high fugacities [15]. The model presents, as we will show, a strong tendency to assume a NaCl-like ordered state at low temperatures. It would be interesting to study the model on the face-centered cubic lattice or another structure that frustrates antiferromagnetic ordering. MEAN-FIELD ANALYSIS The energy of the lattice restricted primitive model is U = 1 2 ∑ i,j sisj rij , (4) where it is understood that the particles occupy distinct lattice sites. We consider the LRPM with independent variables T (temperature) and ρ ≡ N/V . It is helpful to think of this system as a three-state (i.e., spin-1) antiferromagnetic Ising model with long-range interactions (J ∼ 1/r). On a fully occupied lattice (ρ = 1) we characterize ordering by the sublattice charge disparity or staggered magnetization φ. In the context of MFT, this allows us to replace Eq. (4) by the corresponding expression for a nearest-neighbor (NN) lattice gas, multiplied by a suitable factor. Consider first the disordered system. Since the neighbors of any given particle are equally likely to bear the same or opposite charges, the mean-field estimate for the energy is zero, just as for the NN system. Next, suppose there is sublattice ordering; let ρ+ be the fraction of sites in sublattice A occupied by positive particles, etc., and let ρ+ = ρ 2 (1 + φ) ; ρ+ = ρ 2 (1− φ) (5) ρ− = ρ 2 (1− φ) ; ρ− = ρ 2 (1 + φ), (6) so that φ = ρ+ − ρ− ρ , (7) and ρ+ + ρ A − = ρ B + + ρ B − = ρ. (8) Using these, the mean-field estimate for the energy of a pair of nearest-neighbor sites is uNN = ρ A +ρ B + + ρ A −ρ B − − ρ+ρ− − ρ−ρ+ = −ρφ. (9) Similarly, the MF estimate for the energy of a pair of second-neighbor sites (which must belong to the same sublattice) is +ρφ/ √ 2, and so on. Thus MFT gives the energy of the LRPM as UMF = ρφ 2 ∑ i,j (−1)S rij , (10) where S = +1(−1) if sites i and j belong to the same (different) sublattices. The sum (including the factor of 1/2) is simply the electrostatic energy of an ionic crystal. For the simple cubic lattice, 1 2 ∑ i,j (−1)S rij = −3V α 6 , (11) as V → ∞, where α = 1.7457 is the Madelung constant [16]. Thus the meanfield energy per site is uMF = −αρ2φ2/2, while the corresponding expression for a system with nearest-neighbor interactions is −3ρ2φ2. In the context of MFT, then, we may replace the 1/r potential with a nearest-neighbor interaction. From the preceding analysis it appears that we may treat the LRPM, in MFT, as if it were a nearest-neighbor lattice gas, but with an energy (and temperature) scale that is smaller by a factor of α/6 ≃ 0.2913. A more meaningful way of relating the temperature scales, however, is to compare the energy change attending an elementary excitation, in this case, the interchange of a nearest-neighbor pair in a perfectly ordered lattice. In the simple cubic lattice this raises the energy by ∆UNN = 20 (the nearest-neighbor interaction J = 1), since ten nearest-neighbor pairs have like charges after the exchange. For the ionic crystal (NaCl structure) the corresponding energy change is 4(α−1) ≃ 2.9903, so that ∆ULRPM ≃ 0.1495∆UNN . Since the mean-field critical temperature for the (fully occupied) nearest-neighbor lattice gas is 6 on the cubic lattice, the corresponding result for the LRPM is about 0.9. The best numerical estimate for the lattice gas is Tc = 4.5115; the corresponding LRPM value is 0.674. (Extrapolation of our simulation results to ρ = 1 yields Tc ≈ 0.6.) We may now develop the MFT of the LRPM by studying the nearest-neighbor lattice gas, bearing in mind the difference in temperature scales explained above. To estimate the entropy as a function of ρ and φ, we note that the number of allowed configurations on a sublattice of V/2 sites is (V/2)!/[N+!N−!Nv!] where N+ (N−) is the number of positive (negative) particles and Nv = [(V/2) − N+ − N−] the number of vacant sites on the sublattice. Using N± = (1 ± φ)ρV/4, etc., and Stirling’s formula, we obtain the entropy per site, s = −ρ ln ρ− (1−ρ) ln(1−ρ) + ρ ln 2− ρ 2 [(1 + φ) ln(1 + φ) + (1− φ) ln(1− φ)] . (12) Let f = u− Ts be the Helmholtz free energy per site. To investigate the possibility of a free energy minimum with nonzero sublattice ordering φ, we note that